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ABSTRACT 

This paper discusses the identification of fractional- and integer-order systems using the concept of 
continuous order-distribution. Based on the ability to define systems using continuous order-distributions, 
it is shown that frequency domain system identification can be performed using least squares techniques 
after discretizing the order-distribution. 


INTRODUCTION 


Fractional Order Systems 

Fractional order systems, or systems containing fractional derivatives and integrals, have been studied by 
many in the engineering area (Heaviside, 1922; Bush, 1929; Goldman, 1949; Holbrook, 1966; Starkey, 
1954; Carslaw and Jeager, 1948; Scott, 1955; and Mikusinski, 1959). Additionally, very readable 
discussions, devoted specifically to the subject, are presented by Oldham and Spanier (1974) and Miller 
and Ross (1993) and Pudlubny (1999). It should be noted that there are a growing number of physical 
systems whose behavior can be compactly described using fractional system theory. Of specific interest to 
electrical engineers are long lines (Heaviside, 1922), electrochemical processes (Ichise, Nagayanagi, and 
Kojima, 1971; Sun, Onaral, and Tsao, 1984), dielectric polarization (Sun, Abdel wahab, and Onaral, 1984), 
colored noise (Mandelbrot, 1967), viscoelastic materials (Bagley and Calico, 1991; Koeller, 1984; Koeller, 
1986; Skaar, Michel, and Miller, 1988), and chaos (Hartley, Lorenzo, and Qammar, 1995). 

For unknown systems, system identification has become the standard tool of the control engineer. 
Identifying a given system from data becomes more difficult, however, when fractional orders are allowed. 
For integer order systems, once the maximum order of the system to be identified is chosen, the parameters 
of the model can be optimized directly. For fractional order systems, identification requires the choice of 
the number of fractional operators, the fractional power of the operators, and finally the coefficients of the 
operators. Thus, the loss of integer order has significantly complicated the identification process. Previous 
work in this area has been limited (Tsao, Onaral, & Sun, 1989), (Sun Onaral, & Tsao, 1984). These authors 
identify electrode-electrolyte polarization behavior using frequency domain techniques for specifically 
chosen transfer function forms. In what follows, a more general identification technique is presented that 
directly determines the form for the transfer function, any fractional or integer order terms, and the 
coefficients of the individual terms. 
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For reference, we will give the uninitialized Riemann-Liouville definition of the fractional integral, 


0 d;‘>f(t) = -±-\(t-T)‘>- i f(T)dT q>0. (1) 

T(q)l 

and then the fractional derivative is defined as the integer derivative of a fractional integral, 

o = P >o. m 

T(\-p) dt J 

In what follows, it will be important to use the Laplace transform of the fractional integrals and derivatives. 
They are given below, where it is assumed that the initialization is zero, 

L{ 0 d? f(t)\= s* F(s) for all q . (3) 


Generalized Fractional Order Systems 

Here we will use a mechanical system to introduce the idea of continuous order-distributions. The idea will 
be developed slowly, by adding terms individually, and then inductively jumping to continuous order- 
distributions. The usual spring-mass-damper mechanical system is used as a familiar starting place, and is 


given below 

m 0 d;x{t) + b o djx(t) + kx(t) = f(t) (4) 

where x(t) is position of the mass, f{t) is the forcing function on the mass, m is the mass, b is the 
damping, and k is the restoring force. In Laplace transform notation, this can be written as 

(ms 2 + bs + A') X(s) = F(s) . (5) 

Now, it is well known that viscoelastic elements yield fractional order behavior over a wide range of 
frequencies (Bagley and Calico, 1991). Such a viscoelastic element would be described by 

k q 0 d?x(t) = f(t) (6) 

where 0 < q < 1 . Its Laplace transformed representation is 

k q s«X(s) = F(s). (7) 

Now adding this viscoelastic element into the original mechanical system causes the system to become 

(ms 2 + bs + k q s q + A:) X(s) = F(s) . (8) 


The analysis of systems such as these has been discussed at length elsewhere (Lorenzo and Hartley, 1998). 
It is known that the order of the fractional viscoelastic element can truly be anywhere between zero and 
one. Thus it poses no complication to add another viscoelastic element to our system whose order is 
different from the first viscoelastic element. Doing so, the system representation f is 

(ms 2 + bs + k q2 s g '- + + k) X(s) = F(s) . (9) 
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Now it is further known that there are viscoinertial systems that behave as systems with order between one 
and two. Adding two of these elements to our system yields the Laplace domain representation as 

(ms 2 + k (h s‘ u + k qj s‘ h +bs + k q j“- + + k) X(s) = F(s). (10) 


It should be clear that this process could be continued, so that one could express the system equation as 


( n 


]x(s) = F(s) 


J 


(ID 


where 0 < q n <2 for a mechanical system, and N is any positive integer. 


Some materials can display complex thermorheological behavior (Bagley, 1991). This means that the 
order of the viscoelastic element, q n , depends upon the temperature of the material. A large sample of this 
material can be considered as a collection of layered individual elements. If this material is subject to a 
temperature distribution, a corresponding order-distribution will exist throughout the material, as each 
individual element will have its own order. In the limit as the elements get smaller, this leads to the new 
concept of continuous order-distributions. Assuming that we can obtain a material whose order can change 
from zero to two, in the general mechanical system of Equation (11), the summation can be replaced by an 
integral over the system order, 


r 2 'N 

j^)* 9 dq 

vO 


X(s) = F(s) . 


( 12 ) 


While this idea has been introduced using a material example, it is otherwise conceptually simple to jump 
from the summation in Equation (11) to the integral in Equation (12). Thus we now have a very general 
formulation for representing dynamic systems. 

To demonstrate that familiar equations can be written in this form, the common second order mechanical 
system of Equation (5) can also be written using Equation (12), 


( 2 


j[m8(q- 2) + bd(q-l) + k8(q)]s (I .dq 


X(s) = F(s ). 


(13) 


Allowing the restriction on the maximum possible system order to be relaxed from second order, the 
general system representation becomes 


]«*> 


s q dq 


X(s) = F(s). 


(14a) 


It should be observed that this may be thought of as a continuum power series. It also has the time domain 
representation 


J k(q) o d‘'x(t) dq = f(t) . (14b) 

o 
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Clearly, a system could also be represented by a continuum asymptotic series, 


oo 

J*(g)s' 9 dq 


X(s) = F(s). 


( 15 ) 


Equations (14a) and (15) can be combined to give 

A 


dq 


X(s) = F(s). 


(16) 


Also, a more general form, resembling a standard transfer function, is 

foo \ /'oo A 


jk(q)s q dq X(s) = ljg(p) s p dp 


F(s) . 


(17) 


The representations in these last three equations are not considered further in this paper. Equation (14a) will 
be considered at length in the following sections. The next section considers some specific systems that can 
be represented by Equation (14a). After that, the following section develops a system identification method 
based on Equation (14a). 


Analysis of Some Systems with Continuous Order-Distributions 

It is interesting to reconsider the system of Equation (14), by rewriting the exponential in s . That is 


\k(q)e^ s) 

, 0 


A 


dq 


X(s) = P(s)X(s) = F(s). 


) 


(18) 


This equation is effectively a Laplace transform of the function k(q) with the new Laplace variable being 
r = -ln(s). As long as the order-distribution, k(q) , is of exponential order, then the resulting transfer 
function, P(s ), is easy to calculate using this Laplace transform. This calculation is done for several 
mechanical systems in Table 1 and for some other order distributions in Table 2. 


Table 1 presents the transfer functions, P(s ), for some second order systems with continuous order- 
distributions. The first column contains a plot of the order-distributions, k(q)versus q . The second column 
contains the Laplace r -transformed representation using Equation (18). The third column gives the 
corresponding Laplace transfer function for these systems obtained by replacing everywhere in the second 
column r = -ln(s) . 


Table 2 is similar to Table 1, however it contains order-distributions that vary continuously from zero 
to infinity. 

It should be pointed out that systems analysis with continuous order-distributions is somewhat difficult at 
this time. The frequency domain approach is probably the easiest and most reliable, as it is simple to 
replace 5 everywhere by jco , and then vary the frequency. Some frequency responses are given in the 
discussion of the next section. 
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Identification of Order-Distributions from Frequency Domain Data 

Referring back to the general system of Equations (14a) and (18), 


r oo 

dq 


X(5) = F(s) 


it is necessary to first get the transfer function 


X(s) 

F(s) 


1 

OO 

jk(q)s q dq 

o 


1 

P(s) 


= G(s). 


Inverting this equation gives 


\k{q)s q dq = P{s) 
0 


1 

G(s) 


( 19 ) 


( 20 ) 


( 21 ) 


Now, let it be assumed that we have an unknown system, and that a measured frequency response, G(jco ) , 
is available from it. We will assume that the system can be expressed by Equation (14a). That being the 
case, in Equation (21), we can replace s everywhere by jco . Doing so yields, 


\kmm q dq = -^ 

J G( n 


(22) 


The identification problem becomes one of determining the order-distribution, k(q ) , given G(ja >) . As an 
analytical approach is not obvious, a numerical solution that approximates the integral is taken. Many 
approximations for integrals are available, but a simple approach is taken for now and is given below. 

Assuming that the order-distribution, k(q ) , eventually goes to zero as q gets bigger, the integral 
converges and can be replaced by an Euler approximation, or right looking rectangles. The integral of 
Equation (22) then becomes 


s 


t, (jio)" Q Q = 


1 

G(jco) 


(23) 


where Q is the constant sample width in the variable q (a non-constant sample width could also be used), 
and k n is the height of the sampled order-distribution. Remembering that we usually have sampled 
frequency response data. Equation (23) must be satisfied for each data point. Thus we can write 




(jcOj) nQ Q = 


1 

G(jcoj) 


for any j . 


(24) 
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This equation can now be written out for each frequency point 


Qko+Qki ( ;ft), f +Qk 2 (j ft), ) 2Q +Qk 3 {jco t ) 3Q +---+Qk N (;o), f' Q = 

G(jco x ) 

Qk 0 + Qk , Oft), f + Qk, ( /ft) , ) 2(? +fi* 3 Oft), ) 3e +---+0*, v O'ffli = ' (25) 

G(;ft) 2 ) 


Qk 0 + Qk i (ja> M ') +Qk 2 (j(0 M ) +Qk 2 {j(0 M ') ■+ \-Qk N (jco M ) — 


These equations can now be written in matrix form as 


i (M) c 

(M) 22 • 

•• (M) nq ' 

V 


" 1 / G(yft), ) 

l Ua 2 ) Q 

O'a), ) 2Q • 

- Uco 2 ) nq 

k } 

= 

1 / G(jco 2 ) 

.1 

(M/ ) 2Q • 


_k N _ 


l/G(ja„)_ 


(26) 


where it is assumed that the number of frequency response samples, M , is greater than or equal to the 
number of order-distribution samples, N . This equation can be written more compactly as 

Wk = g. (27) 

where W is the big matrix including Q , and k and g are the vectors in Equation (26). Clearly, if 

M > N , then a least squares solution, or matrix pseudoinverse, can be used to solve for the sampled order- 
distribution, k(q) . 

One problem with this approach is that the matrix W tends to become singular as the order-distribution 
sample size, Q , gets small, and also as the number of order samples, N , gets large. The user might also be 
tempted to use the form of Equation (18) which more closely resembles a Fourier transform in log- 
frequency. Although solvable, most software will have some difficulties with the terms 

nQln(jO) m ) 

as it is hard to stay on the primary Riemann sheet, after multiplying the log term by a number that can 
become large. This was not done in this study. In our studies, we used the Matlab pseudoinverse 
command pinv on Equation (27), as it was numerically more reliable than performing the pseudoinverse 

using ^ = (w rr VV r ) W T g . Pre-multiplication on both sides by a scaling matrix was also performed, and 
that also improved the conditioning of the matrix. 
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EXAMPLE APPLICATIONS 


Several examples are now presented to demonstrate the utility of the identification process. The 5-domain 
transfer functions from the tables are used to calculate a numerical frequency response by replacing 
5 = jco . These frequency responses are then used to generate the g vector of Equation (26). 


In every case but one, Equations (26) were used with Q = 0.1 , N = 25, and co logarithmically distributed 

between 10” 2 to 10 2 . The examples chosen are entries 1, 2, 3, 7, and 8 from Table 1, entry 1 in Table 2 
(with the exponent being 2q and Q=0.5, see Figure (6)), as well as the following 


G(s) = 


1 

5 2 +1.45 1 - 5 +5-fl.45°- 5 +l 


(28) 


and 


G(s) = 


1 

5 ~ + 0.55 - 1 - 1 


(29) 


The identified distributions, &(<?) versus <7 , and their frequency magnitude responses, in Bode form, are 
given in Figures 1 through 8. It should be noted in the frequency responses, that both the original and 
identified frequency response plots lie on top of one another, with negligible error, for the systems 
presented. The phase responses were not included, but were similarly accurate. In looking at the identified 
order distributions, an interesting oscillation is observed whenever there is a discontinuity in order, 
somewhat reminiscent of the Gibb’s phenomenon from Fourier analysis. It was found that using different 
data reconstructors (filters associated with the sampling process) in the q -variable it was sometimes 

possible to reduce this. 

The identification process outlined above yields as its result k(q ) , the order-distribution, for an input 
frequency response G(jco). To obtain the original transfer function, G(s ) , from k(q) (the identified order- 
distribution), the following expression must be used 

g(s)b^ , (30) 

t,k n s nQ Q 

n~0 

which is a sampled form of equation (20). This is particularly important for the results pertaining to 
Equations (28) and (29) that have discrete order-distributions. In these cases, it should also be observed that 
the order spikes that one would expect are somewhat smeared in order (see Figures (7) and (8)). 
Techniques to concentrate the k{q) distribution peaks into fewer specific (discrete) q terms need to be 
carefully developed. These issues remain to be resolved. 


DISCUSSION 

This paper develops an identification method based on the concept of continuous order-distributions. This 
technique allows the identification of both standard fractional and integer order systems containing 
individual, or discrete, terms, as well as systems with continuous order-distributions. The technique was 
demonstrated on systems with both continuous order-distributions and discrete order-distributions. The 
effect of noise, numerical truncation, and expected accuracy of the results has not been studied at this time. 
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Table 1 . — Order-distributions, and their transfer functions for some systems up to order 2. 


Order Distribution Laplace r-transform, Eqn. (18) Laplace transfer function, r=-ln(s) 
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Table 2. — Order-distributions and their transfer functions for several systems with 
order continuously distributed up to infinity. 


Order Distribution for Material 
k(q) = e~ q 


k{q) = — 
q + a 


k(q) = ^(-l) n+] 8(q-n) 

n= 1 


k(q)=sm 2 (aq) 


k(q) = -y-\og(q) 


Laplace r-transform, Eqn. (18) Laplace transfer function, r=-ln(s) 


P(r) = 


r+1 


P(s) = 


1 

l-ln(j) 


( 1 ) 


P(r)--e ar Ei(-ar) 


P(s) = -aEi(-a ln(s)) (2) 


P(r) = 


1 

si nil (r) 


P(s) = 


1 

sinh(ln(s)) 


(3) 


P(r) = 


la 2 

r(r 2 +4a 2 ) 


P(s) = 


— la 


ln(^) (ln(s)) 2 +4a 2 


(4) 


P(r) = 


ln(r) 

r 


P(s) = 


ln(-ln(s)) 

-ln(s) 


(5) 
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Figure 1. — Identified order-distribution and magnitude frequency response for entry 1 in Table 1. 



q co rads 


Figure 2. — Identified order-distribution and magnitude frequency response for entry 2 in Table 1. 



Figure 3. — Identified order-distribution and magnitude frequency response for entry 3 in Table 1. 


NAS A/TM — 1 999-209640 


12 



















Figure 7. — Identified order-distribution and magnitude frequency response for text Equation (28). 
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Figure 8. — Identified order-distribution and magnitude frequency response for text Equation (29). 
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